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Abstract 

A method for simulating power law noise in clocks and oscillators is 
presented based on modification of the spectrum of white phase noise, 
then Fourier transforming to the time domain. Symmetric real matrices 
are introduced whose traces-the sums of their eigenvalues-are equal to the 
Allan variances, in overlapping or non-overlapping forms, as well as for 
the corresponding forms of the modified Allan variance. Diagonalization 
of these matrices leads to expressions for the probability distributions for 
observing a variance at an arbitrary value of the sampling or averaging 
interval r, and hence for estimating confidence in the measurements. A 
number of applications are presented for the common power-law noises. 

1 Introduction 

The characterization of clock performance by means of average measures such as 
Allan variance, Hadamard variance, Theo, and modified forms of such variances, 
is widely applied within the time and frequency community as well as by most 
clock and oscillator fabricators. Such variances are measured by comparing 
the times tk on a device under test, with the times at regular intervals fcro on a 
perfect reference, or at least on a better reference. Imperfections in performance 
of the clock under test are studied by analyzing noise in the time deviation 
sequence x k = tk— kro, or the fractional frequency difference during the sampling 
interval r = stq: 

A M = ( x k+s ~ X k )/ (ST ). (1) 

The frequency spectrum of fractional frequency differences can usually be ad- 
equately characterized by linear superposition of a small set of types of power 
law noise. The frequency spectrum of the fractional frequency differences of a 
particular noise type is given by a one-sided spectral density 

S y (f) = h a f a , />0. (2) 
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(The units of S y (f) are Hz -1 .) For the common power-law noise types, a 
varies in integral steps from +2 down to -2 corresponding respectively to white 
phase modulation, flicker phase modulation, white frequency modulation, flicker 
frequency modulation, and random walk of frequency. 

Simulation of clock noise can be extremely useful in testing software algo- 
rithms that use various windowing functions and Fourier transform algorithms 
to extract spectral density and stability information from measured time devi- 
ations, and especially in predicting the probability for observing a particular 
value of some clock stability variance. This paper develops a simple simulation 
method for a time difference sequence that guarantees the spectral density will 
have some chosen average power law dependence. Expressions for the common 
variances and their modified forms are derived here that agree with expres- 
sions found in the literature, with some exceptions. This approach also leads 
to predictions of probabilities for observing a variance of a particular type at 
particular values of the sampling time. A broad class of probability functions 
naturally arises. These only rarely correspond to chi-squared distributions. 

This paper is organized as follows. Sect. 2 introduces the basic simulation 
method, and Sect 3 applies the method to the overlapping Allan variance. Sect. 
4 shows how diagonalization of the averaged squared second-difference operator, 
applied to the simulated time series, leads to expressions for the probability of 
observing a value of the variance for some chosen value of the sampling or 
averaging time. Expressions for the mean squared deviation of the mean of 
the variance itself are derived in Sect. 6. The approach is used to discuss the 
modified Allan variance in Sect. 7, and the non-overlapping form of the Allan 
variance is treated in Sect. 8. Appendix 1 discusses evaluation of a contour 
integral for the probability, first introduced in Sect. 4, in the general case. 

2 Discrete Time Series 

We imagine the noise amplitudes at Fourier frequencies f m are generated by a 
set of N normally distributed random complex numbers w n having mean zero 
and variance a, that would by themselves generate a simulated spectrum for 
white phase noise. These random numbers are divided by a function of the 
frequency, |/ m | A , producing a spectral density that has the desired frequency 
characteristics. For ordinary power law noise, the exponent A is a multiple 
of 1/2, but it could be anything. The frequency noise is then transformed 
to the time domain, producing a time series with the statistical properties of 
the selected power law noise. The Allan variance, Modified Allan variance, 
Hadamard Variance, variances with dead time, and other quantities of interest 
can be calculated using either the frequency noise or the time series. 

In the present paper we discuss applications to calculation of various versions 
of the Allan variance. Of considerable interest are results for the probability 
of observing a value of the Allan variance for particular values of the sampling 
time t and time series length N. The derivations in this paper are theoretical 
predictions. A natural frequency cutoff occurs at fh = 1/(2t ), where t is 
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the time between successive time deviations. This number is not necessarily 
related in an obvious way to some hardware bandwidth. The measurements are 
assumed to be made at the times fcro, and the time errors or residuals relative 
to the reference clock are denoted by Xk- The averaging or sampling time is 
denoted by r = st , where s is an integer. The total length of time of the 
entire measurement series is T = Nt . The possible frequencies that occur in 
the Fourier transform of the time residuals are 

Noise Sequences. In order that a set of noise amplitudes in the frequency 
domain represent a real series in the time domain, the amplitudes must satisfy 
the reality condition 

w- m = (w m )*. (4) 

N random numbers are placed in N/2 real and N/2 imaginary parts of the 
positive and negative frequency spectrum. Thus if w m = u m + iv m where u m 
and v m are independent uncorrelated random numbers, then (w m )* = u m — iv m . 
Since the frequencies ±l/2r represent essentially the same contribution, v N / 2 
will not appear. We shall assume the variance of the noise amplitudes is such 
that 

{(w m )*w n ) = {u 2 + v 2 )S mn = 2a 2 S mn ; m^0,N/2. (5) 

Also, < >=< u 2 — v 2 + 2iuv >= for m ^ N/2. The index m runs from 
—N/2 + 1 to N/2. In order to avoid division by zero, we shall always assume 
that the Fourier amplitude corresponding to zero frequency vanishes. This only 
means that the average of the time residuals in the time series will be zero, and 
has no effect on any variance that involves time differences. 

We perform a discrete Fourier transform of the frequency noise and obtain 
the amplitude of the k th member of the time series for white PM: 

N/2 

x k = -t= 2^ e N Wm - w 

* N m=-N/2+l 

The factor Tq is inserted so that the time series will have the physical dimensions 
of time if w m has the dimensions of frequency. We then multiply each frequency 
component by |/o// m | A - This will generate the desired power-law form of the 
spectral density. The time series will be represented by 

2 N / 2 I f |A 

V N m= _JV/2+l l/m| 

The constant factor |/o| A has been inserted to maintain the physical units of 
the time series. The noise level is determined by f . For this to correspond to 
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commonly used expressions for the one-sided spectral density, Eq. @, we shall 
assume that 

T 2 |/0| A _ / K (8) 



N V IQ-k 2 <t 2 {Nt ) 

We shall show that if 2 A = 2 — a the correct average spectral density is obtained. 
The simulated time series is 



2Tvimk 

e N 



Xfc = Vi6^V^)?17^- (9) 

The average (two-sided) spectral density of the time residuals is obtained from 
a single term in Eq. 

167rV 2 (ATTo)/^\ A/ / ~ 8^ 2 / 2 



where A/ = 1/(JVto) is the spacing between successive allowed frequencies. The 
average (two-sided) spectral density of fractional frequency fluctuations is given 
by the well-known relation 

s y (f) = (27rf) 2 s x (f), (11) 
and the one-sided spectral density is 

s yW = {°2s y (f) = h a r, />o, (12) 

where 2A = 2 — a. 



3 Overlapping Allan Variance 

Consider the second-difference operator defined by 

V2r 2 

The fully overlapping Allan variance is formed by averaging the square of this 
quantity over all possible values of j from 1 to N — 2s. Thus 

/ 1 N - 2s \ 

In terms of the time series, Eq. (|9|), the second difference can be reduced using 
elementary trigonometric identities: 
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A (2) = 



327r 2 r 2 CT 2 (iVTo) 

27riTn(j + 2s) 2 J iTTi(j + a) 2.ii.(j) 



^(e- . -2,- . + e 



27r 2 rV 2 (iVTo) ^ |/ m P V # 

( 2) 

We form the averaged square of A) j by multiplying the real quantity times its 
complex conjugate, then averaging over all possible values of j 



ay{T) ~ 2TT 2 r 2 (NTn)(N -2s) L>\ a 



sin (^T) sin 
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2tt^(Nt )(N-2s) a 2 / |/ ro / n |* 

2ir»(m-n)(j + a) 

x e « . (16) 

The average of the product of random variables only contributes 2cr 2 when 
m = n (see Eq. (JSJ) , except when m = n = N/2 where (\w N / 2 \ 2 ) = v 2 ■ The 
Allan variance reduces to 

- f y v :^)) < i i ( si °<«)\ 

» U ^(JVro)^ | /m px +2 (/jv/2) « J- 

since every term in the sum over j contributes the same amount. The zero 
frequency term is excluded from the sum. For convenience we introduce the 
abbreviation 

A ' = ;sw (18) 

If we sum over positive frequencies only, a factor of 2 comes in except for the 
most positive frequency and so 

4 



»/2-l (sinpjfi)) Ysinf) 



* , = n,S^r^ + J wkY (19) 

If the frequencies are spaced densely enough to pass from the sum to an integral, 
then A/ = (iVro)" 1 and 

^-E F (i/™i)^ / ( 2 °) 
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and we obtain the well-known result [6 

-»-^J Q fh ^f(^fr))\ (21) 

Similar arguments lead to known expressions for the non-overlapping version 
of the Allan variance as well as for the modified Allan variance. These will be 
discussed in later sections. 



4 Confidence Estimates 

In the present section we shall develop expressions for the probability of ob- 
serving a particular value A Q for the overlapping Allan variance in a single 
measurement, or in a single simulation run. A a is a random variable repre- 
senting a possible value of the overlapping variance. We use a subscript "o" 
to denote the completely overlapping case. To save writing, we introduce the 
following abbreviations: 

|/m|* \ N ) 

ni _ ( sin (V)) ( 2nm(j + s) 
Gm ~ U™F &m l N 

The dependence on s is suppressed, but is to be understood. We write the 
second difference in terms of a sum over positive frequencies only, keeping in 
mind that the most positive and the most negative frequencies only contribute 
a single term since sin(7r(j + s)) =0. The imaginary contributions cancel, and 
we obtain 

A ^ = ^E(^v + G ™v)- (23) 

m>0 

There is no term in %^ 2 - It is easy to see that the overlapping Allan variance 
is given by 

°» = E E ((^) 2 + ( G ^) 2 ) • (24) 

To compute the probability that a particular value A a is observed for the Allan 
variance, given all the possible values that the random variables Ux, Hi, ...Ujj/2 
can have, we form the integral 

™ - 1 { A ° -f^m> ! ) n ( e -^^) . (25, 

J \ j ' m >0 v ' 

The delta function constrains the averaged second difference to the specific 
value A a while the random variables u\, v\, ...u m , v m , ...u^/2 range over their 




G 



(normally distributed) values. There is no integral for vm/2- Inspecting this 
probability and the Eq. (|23|) for the second difference indicates that we can 
dispense with the factors of <r -1 and work with normally distributed random 
variables having variance unity. Henceforth we set a = 1. 

The exponent involving the random variables is a quadratic form that can 
be written in matrix form by introducing the N — 1 dimension column vector 
U (the zero frequency component is excluded) 

U T = [uivi ...u m v m ,... wat/2- i,u n/2 ] . (26) 

Then 

\Y, ( ^n + <) = \u T U=\u T lU, (27) 

m>0 

where 1 represents the unit matrix. The delta-function in Eq. (|25p can be 
written in exponential form by introducing one of its well-known representations, 
an integral over all angular frequencies w:[S] 

P(A o) = I" *Le*>(*>-*±X* W) J] (e-^*g*¥) . (28) 

The contour of integration goes along the real axis in the complex u> plane. 

The squared second difference is a complicated quadratic form in the random 
variables u\,vi, ...u m ,v m , ...u N / 2 . If this quadratic form could be diagonalized 
without materially changing the other quadratic terms in the exponent, then the 
integrals could be performed in spite of the imaginary factor i in the exponent. 
To accomplish this we introduce a column vector C J that depends on j, m, s, N 
and whose transpose is 

(C 3 ) = [F( , G\, ...F^, G 3 m , ...G^ v y 2 _ 1 , i*jy/ 2 ] • (^^) 

Dependence on s is not written explicitly but is understood. The column vec- 
tor has N — 1 real components. It contains all the dependence of the second 
difference on frequency and on the particular power law noise. We use indices 
{m, n} as matrix (frequency) indices. The (scalar) second difference operator 
can be written very compactly as a matrix product 

= VKi&fU = ^KU T C J . (30) 

Then 

if^£( 4 2)'- or (ir5£' ! '(^) R (31) 

j v 7 x 3 7 

The matrix 

H ° = WTYX Cj{Cj)T (32) 

j 

is real and symmetric. H Q is also Hcrmitian and therefore has real eigenvalues. A 
real symmetric matrix can be diagonalized by an orthogonal transformational, 
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[2] which we denote by O. Although we shall not need to determine this orthog- 
onal transformation explicitly, it could be found by first finding the eigenvalues 
e and eigenvectors tp e of H Q , by solving the equation 

H a ^ = eA . (33) 

The transformation O is a matrix of dimension (JV — 1) x (N — 1) consisting of 
the components of the normalized eigenvectors placed in columns. Then 

H a O = OE , (34) 

where E is a diagonal matrix with entries equal to the eigenvalues of the matrix 
H Q . Then since the transpose of an orthogonal matrix is the inverse of the 
matrix, 

T H a O = E . (35) 

The matrix H Q is thus diagonalized, at the cost of introducing a linear trans- 
formation of the random variables: 

^ E ( A 2) 2 = UTH ° U = U t OO t H OO t U = (U T 0)E(0 T U) . (36) 

3 

We introduce N — 1 new random variables by means of the transformation: 

V = T U. (37) 
Then the term in the exponent representing the Gaussian distributions is 

1 1 1 ^ N ~ 1 

- -U T 1U = --U T OW T U = --V T 1V = -- V V? . (38) 

2 2 2 2 ^ " v ; 

n— 1 

The Gaussian distributions remain basically unchanged. 

Further, the determinant of an orthogonal matrix is ±1, because the trans- 
pose of the matrix is also the inverse: 

detiO^O) = 1 = det(0 T 0) = (det(O)) 2 . (39) 

Therefore, changes in the volume element are simple since the volume element 
for the new variables is 



dV 1 dV 2 ...dV N - 1 



dot 



dV„. 



dU 1 dU 2 ...dU N - 1 



dU n , 

= \det(0)\dUidU 2 -dU N - 1 

= dU 1 dU 2 ...dU N - 1 . (40) 

After completing the diagonalization, 

^E(^) 2 -E^ 2 - (4i) 
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The probability is therefore 



P(A o) = J %H A °^ n (e-* • (42) 

An eigenvalue of zero will not contribute in any way to this probability since 
the random variable corresponding to a zero eigenvalue just integrates out. 

Let the eigenvalue have multiplicity \Xi , by which is meant that the eigen- 
value ej is repeated times. Integration over the random variables then gives 
a useful form for the probability: 

/+oo j p iuA 
_ s n, ( i +aw ,)>/. - <«» 

Finally the contour integral may be deformed and closed in the upper half 
complex plane where it encloses the singularities of the integrand. This is dis- 
cussed in Appendix 1. Knowing the probability, one may integrate with respect 
to the variance to find the cumulative distribution, and then find the limits 
corresponding to a 50% probability of observing the variance. 

Properties of the eigenvalues. First, it is easily checked that the probability 
is correctly normalized by integrating over all A a and using properties of the 
delta-function: 



I Pi A )dA f °° d " I^±**> f °° 

J { ° ] ° loo 27T FL(1 + 2*6^/2 - J_ 



+ 2teiw)w/2 

/+oo 
du>5{ui) = 1 . (44) 
-co 



Second, the eigenvalues are all either positive or zero. The eigenvalue equation 
for the eigenvector labeled by e is: 

_JL_Y,Ci{C>) T ^=e^. (45) 

3 

Multiply on the left by ipj; assuming the vector has been normalized, we obtain 

2 



K 



N -2s 

j 



£(W^) >o. 



(46) 



Thus every eigenvalue must be positive or zero. 

Next let us calculate the trace of H Q . Since the trace is not changed by an 
orthogonal transformation, 

Trace(0 T ff 0) = Tracc(ff 00 T ) = Tracc(ff o 0C>- 1 ) 

= Tracc(tf ) = ^ e, . (47) 
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The sum of the diagonal elements of H a equals the sum of the eigenvalues of 
H a . If we then explicitly evaluate the sum of the diagonal elements of H a we 
find 

* 3 

K EE ( (F^) 2 + (Ginf) 



N-2s 

j m>0 



sin 2$2 



. n \Jm\ 
m>0 

Every term labeled by j contributes the same amount. We obtain the useful 
result that the overlapping Allan variance is equal to the sum of the eigenvalues 
of the matrix H a . Similar results can be established for many of the other types 
of variances. 

Distribution of eigenvalues. The eigenvalue equation i?oVv — etp^ produces 
many zero eigenvalues, especially when r is large. The dimension of the matrix 
H a is therefore much larger than necessary. Numerical calculation indicates 
that for the completely overlapping Allan variance, the eigenvalue equation has 
a total of N — 1 eigenvalues, but only N — 2s non-zero eigenvalues; the number 
of significant eigenvalues is in fact equal to the number of terms in the sum over 
j in the equations: 

N-2s 

^ E E (c m y(c n y^ n = ^ m . (49) 

n j 

The factorized form of H Q , that arises on squaring a difference operator, permits 
the reduction of the size of the matrix that is to be diagonalizcd. We introduce 
the quantities 

<Ai = £( c ") j ^- (50) 

n 

We are using the Greek index fi to label a non-zero eigenvalue and the index v 
to label a zero eigenvalue. The eigenvalue equation becomes 

^_^(C m )^ = e^. (51) 

3 

Multiply by (C m ) 1 and sum over the frequency index m. Then 

-E—J2(C m ) l (C m y^ = e^. (52) 

m,j 

This is an eigenvalue equation with reduced dimension N — 2s rather than N—l, 
since the number of possible values of j is N — 2s. The eigenvalue equation can 
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be written in terms of a reduced matrix H re d, given by 

m 

The indices i, j run from 1 to N — 2s. Eigenvalues generated by Eq. (|52| are all 
non-zero. To prove this, multiply Eq. (|52j) by (j> 1 and sum over We obtain 



N -2s 

m v I 



The eigenvalue cannot be zero unless 

£(C m )'^ = (55) 
I 

for every m. The number of such conditions however is larger than the number 
N — 2s of variables, so the only way this can be satisfied is if <j> 1 = 0, a trivial 
solution. Therefore to obtain normalizable eigenvectors from Eq. (j52|) . the 
corresponding eigenvalues must all be positive. This is true even through some 
of these conditions may be trivially satisfied if the factor sm.{wms / N) vanishes, 
which happens sometimes when 

ms = MN (56) 

where M is an integer. Every time a solution of Eq. (|56[) occurs, two equations 
relating components of <^ are lost. Suppose there were n solutions to Eq. (|56[) : 
then the number of conditions lost would be 2n. The number of variables is 
N — 2s and the number of conditions left in Eq. (|B3)l would be TV — 1 — 2n. The 
excess of conditions over variables is thus 

N - 1 - 2n- (N- 2s) = 2(s - n) - 1 . (57) 

In Appendix 2 we prove that under all circumstances 2 (a — n) — 1 > 0. 

We temporarily drop the subscript o since the remainder of the results in 
this section are valid for any of the variances. If the eigenvalues are found and 
the appropriate matrix is diagonalized, we may compute the probability for 
observing a value of the overlapping variance, denoted by the random variable 
A by 



dui i,J a-v t f.v\ rr/e V ' 2/2 dV,, 



P {A ) = j°° ±L e ^-v T Ev) JJ 



V2tt 

2^n( 1 + 2 ^ w ) At,/2 



(58) 



Case of a single eigenvalue. If a single eigenvalue occurs once only, the general 
probability expression, Eq. (|43p. has a single factor in the denominator: 

P(A) = ± / . (59) 

2tt J —oo V 1 + 2iwe 
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The integral is performed by closing the contour in the upper half complex u) 
plane. There is a branch point on the imaginary axis at Li = i/(2e), and a 
branch line from that point to infinity. Evaluation of the integral gives: 

1 p -A/(2* 2 y (T)) 

P(A)= . -=— (60) 



This is a chi-squared distribution with exactly one degree of freedom. The 
computation of the confidence interval for a given s is simple. The cumulative 
probability obtained from Eq. ([60]) is 



fA i -x/(2a 2 y (r)) a 

HA)= I T ^ = erf( ). (61) 



The ±25% limits on the probability of observing a value A are then found to 
be 1.323cjy and 0.1015(7^, respectively. An example of this is plotted in Figure 
2. (A is a variance, not a deviation.) 

Case of two distinct non-zero eigenvalues. For the overlapping variance, 
when s has its maximum value N/2 — 1 there are two unequal eigenvalues. The 
probability integral can be performed by closing the contour in the upper half 
plane and gives the expression 

P(A) = -^e^(^)lo(^(I-i)). (62) 

where Iq is the modified Bcssel function of order zero. The probability is cor- 
rectly normalized. It differs from a chi-squared distribution in that the density 
does not have a singularity at A — 0. This is illustrated in Figure 1. 

Evaluation of the contour integral when there are more than two distinct 
eigenvalues is discussed in Appendix 1. If an eigenvalue occurs an even number 
2n times, the corresponding singularity becomes a pole of order n and a chi- 
squared probability distribution may result; this has only been observed to occur 
for white PM. 



5 Variance of the Variance 

One of the goals of this investigation is to estimate the uncertainty in the vari- 
ance when some value of the variance is measured. This can be rigorously 
defined if the probability P(A Q ) is available, but it may be difficult to evaluate 
the integral in Eq. (|4*5|) . In this case, one might be interested in a measure such 
as the rms value of the variance, measured with respect to the mean variance 
<jy (t) . One method for obtaining this measure is obtained from the probability 
without introducing a Fourier integral representation for the delta function: 

A 2 P(A )dA = f A 2 S(A -$>^ 2 )I1 (e- v ^^£\dA 
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/ ( e ^) V* - e n ( e "^ /2 i=) ^ 



e-^/ 2 ^ I . (63) 



Expanding the square in the integrand, there are N fourth-order terms and 
N(N — 1) cross terms so we get 

J AlP(A )dA = Y.^( V ')+ 2 Y.^(V 2 ) 2 

= 3 E e2 + 2 E e ^- ( 64 ) 

To obtain the variance of the variance, we must subtract 

(A ) 2 = ^+ 2 Y,^- ( 65 ) 

The result is 

(^)-(A ) 2 = 2^ £4 2 . (66) 

i 

Therefore the rms deviation of the variance from the mean variance equals y/2 
times the square root of the sum of the squared eigenvalues. This result is not 
very useful if there are only a small number of eigenvalues, since the factor 
y2 can make the rms deviation from the mean of the variance larger than the 
variance. This is because the probability distribution of the variance is not a 
normal distribution, and can lead to unreasonably large estimates of confidence 
intervals. 

If the eigenvalues or the probabilities are not readily available, a similar 
confidence estimate can be obtained from an alternative form of Eq. (|66p by 
considering the trace of the square of H a : 

2 E e2 = 2 E(E( F ™ F ™ + G ™ G ™)) 2 - ( 6? ) 

i i.j m>0 

To prove this, consider diagonalization of H Q by the orthogonal transformation 
O, and compute 

2^e- = 2Trace(£ 2 ) = 2Trace(0 T H o 00 7 H a O) = 2Traceff 2 

i 

= W^r^ E ( E ((c m T(c m y) £ ((c n y(c n y) 

9K 2 f 

= w-^s? £ ( ^ {{c m )\c m y) £ ((c n y(c n y) 

2h ~ E(E(^ F ™ + G - G ™)) ■ ( 68 ) 



(N- 2s) 2 
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If there are not too many terms in the sums over i and j then the sum over the 
frequency index in Eq. (|67|l becomes 

m m>0 m 

This leads to some useful approximation schemes, but these will not be discussed 
further in this paper. 

Usually the variance of the variance is larger than the range of possible 
values of the variance computed from the ±25% limits obtained from calculated 
probabilities. An example is given in Figure 1, where the confidence intervals are 
plotted for the non-overlapping variance for N = 1024. For s ranging between 
342 and 511 there is only one eigenvalue, so the rms deviation of the variance 
from the variance is y/2 times the variance. This gives an upper limit on the 
confidence interval that is 1.414 times the variance, larger than that obtained 
from the actual probability distribution. The medium heavy lines in Figure 1 
show the true ±25% probability limits obtained from calculations such as Eq. 

una 



6 Modified Allan Variance 

The modified Allan variance is defined so that averages over time are performed 
before squaring and averaging. We use a subscript m to distinguish this form 
from the overlapping form of the variance. The definition is 



3=1 

Summing the expression for the second-difference given in Eq. (I15[) . 



l>r^ A (2)_ / ha \ - w m , . irms , 2 



sin^ 



sin — 
bin N 



(71) 



Squaring this, we write the complex conjugate of one factor and obtain the 
ensemble average 



6 

Tims 



2 / \ = ^ v ( Wm <) v 81D n j (79] 

mK} 2sWa*(NT )^ n \f m f n \> (sin^) 2 ' 
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Figure 1: Comparisons of average non-overlapping Allan variance (light line) 
with rms deviations from the mean of the variance (heavy lines) and true ±25% 
limits (medium heavy lines) for N = 1024 data items in the time series. Re- 
sults for the Allan variance for two independent simulation runs are plotted for 
comparison. 
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Noise Type 


Sy(f) 


Mod ^(r) 


A 


a 


White PM 




7r 2 T 3 \± T 1Ss .i) 





2 


Flicker PM 


hif 


3h ± ln(256/27) 
8tt 2 t 2 


1 
2 


1 


White FM 


h 


4r \^ T 2p7 


1 





Flicker FM 


h-if- 1 




3 
2 


-1 


Random Walk 


h- 2 f- 2 




2 


-2 



Table 1: Asymptotic expressions for the Modified Allan Variance, in the limit 
of large sampling times r = stq. 



Using Eq. (JSJand writing the result in terms of a sum over positive frequencies 
only, 



6 

urns 



"' 11 ' sin^i 



Passing to an integral when the frequency spacing is sufficiently small, 

h * 2 S 2 T*P ( Sill (,/ T|) )) 2 

This agrees with previously derived expressions [12]. The appearance of the sine 
function in the denominator of this expression makes the explicit evaluation 
of the integral difficult. In Table 1 we give the leading contributions to the 
modified Allan variance for very large values of the averaging time r. In general 
there are additional oscillatory contributions with small amplitudes that are 
customarily neglected. These results do not agree however with those published 
in [7]; the leading terms do agree with those published in [8]. 

Eigenvalue structure for modified Allan variance. For the modified case, 
instead of Eq. Q32p there are two summations; we use a subscript m on H to 
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denote the modified form: 

^ = 55>'B cfc ) T - (75) 

3 J k 

We then seek solutions of the eigenvalue equation 

j k 

Multiply by ^2[(C l ) T and sum over the frequency index. 

J E ((C m ) l ) T (C m y(J2(C k ) T ^) = e(£(C) T lk) ■ (77) 

l,j,m k I 

The matrix equation has been reduced to a scalar equation for the quantity 

* = £(C*) T &- (78) 

k 

Here the matrix product of (C k ) T with ijj e entails a sum over all frequencies. 
Since we have a scalar eigenvalue equation for e, there can be one and only one 
eigenvalue, which is easily seen to be the same as given by Eq. (|73j) for each 
sampling time r: 

e = f E ((C m )0 T (G m ) fc = ^(r). (79) 

l,k,m 

The probability distribution will be of the form of Eq. (|60p . a chi-squared 
distribution with one degree of freedom. 



7 Non-overlapping Allan variance 

In the non-overlapping form of the Allan variance, the only change with respect 
to Eq. (fT^)) is that the sum over j goes from 1 in steps of s up to j m ax < N — 2s. 
We denote the number of values of j by n max . The average non-overlapping 
variance is the same as that for the overlapping case, but the probability distri- 
butions are different. The matrix H no takes the form 

K Imax 

(H no ) mn = - — E (c m y{c n y. (so) 

max j—i i_^_2s.. 

The subscript "no" labels a non-overlapping variance, and the indices m and n 
label frequencies. Here there are only n max terms in the sum since the values 
of j skip by s. n max is given by: 

N — 1 

n ma x = [ J - 1 , (81) 
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where \Q\ denotes the largest integer less than or equal to Q. In order to 
diagonalize H no and compute probabilities, we look for eigenvalues by seeking 
solutions of: 

^(-ff« )mnV'n = «Am- (82) 
n 

We define a reduced eigenvector by 

0i=£(Cn) J >n. (83) 



The eigenvalue equation reduces to 

K 



3 



^2(c m y^ = «/w (84) 



Multiply this equation by (C m ) 1 and sum over the frequency labels m. Then 
K 



^^((C m y(C m y)^ = e^. (85) 

j m 

Let 

Jij = — ^2(C m ) l (C m y = —(C l ) T Ci; (86) 

m 

the sum over frequencies is accomplished by the matrix multiplication. Just as 
for the overlapping case, if we calculate the trace of the matrix Jij , we find that 
since each term in the sum over j contributes the same amount, 

Trace( Jy ) = £ CjCj = a 2 no (r) . (87) 

'"max 

J 

Because the trace remains unchanged under an orthogonal transformation, the 
non-overlapping Allan variance will be equal to the sum of the eigenvalues of 
Eq. (|85p. The probability functions for a given r can differ from those for 
the overlapping or modified cases because the number of eigenvalues and their 
multiplicities may be different. The eigenvalues are still all greater than or equal 
to zero. 

If the eigenvalues are found and the matrix Jij is diagonalized, we may com- 
pute the probability for observing a value of the overlapping variance, denoted 
by the random variable A nol by 

p{Ano)= L^ e { 

2^n( 1 + 2 *e<w)' 1 '/ 2 ' (88) 



18 



For example, when there are only two distinct non-zero eigenvalues, The matrix 
Jik will have elements 



•h,l — Jl+s.l + s — 7: a no'i 



1 

2' 

•h,i+s = Ji+.,i = I" ^|-(sin(^r/)) 4 cos(27rr/) . (89) 
Jo K^tj) 

The eigenvalues are then 

£2 = \<J 2 no - | Ji,i+.| . (90) 



The probability is of the form of Eq. (1621) . Figure 1 shows an example of this for 
flicker FM noise, for N — 64 items in the time series. The variance corresponding 
to t = 19to was extracted from each of 4000 independent simulation runs and 
a histogram of the resulting values was plotted. Good agreement with the 
predicted probability distribution can be seen. 

8 Other variances 

The analysis methods developed in this paper can be extended to other vari- 
ances. For example, the Theo variance|13j is defined by an overlapping average, 

^ N—m . 
<T 2 Th (s,T Q ,N) = j^— s E ^2 X 

2—1 

s/2-l 1 . y 
^ g/2 - 8 ~ Xi - S + S / 2 ^ ~ ( X i+s ~ X i+8+8/2)j ■ (91) 

where s is assumed to be even, and t = Sstq/A. For any of the power law noises, 
and after passing to an integral, this expression can be transformed with the 
aid of elementary trigonometric identities to 

, , Vl r h Sy(f)df 2T " 

^Th{ T , T 0,^) 
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r s y (f)df '"w . 4 y 

J -J^JyT X. -^ 1 n( 7 r/ K ro)sm( 7 r/(-r-Kro))j 



(92) 

The integer k in the denominator makes this variance more difficult to evaluate, 
but a factorized quadratic form can still be constructed and the eigenvalue 
structure and resulting probabilities can be analyzed. 

The Hadamard variance is defined in terms of a third difference, and is widely 
used to characterize clock stability when clock drift is a significant issue[]3]. A 
third difference operator may be defined as 

Aj> = /g^j ( x j+3s - 3X j+2s + X j+S - Xj) . (93) 
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Figure 2: Simulation of Allan variance with N = 64, s = 19. For this case 
there are two distinct eigenvalues in the matrix for the non-overlapping case. 
The variance for s = 19 was extracted from each of 4000 independent runs and 
a histogram of the values obtained was constructed for comparison with the 
probability, Eq. (I6"2"j) . Chi-squared distributions with 1, 2, and 3 degrees of 
freedom are plotted for comparison. 
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The completely overlapping Hadamard variance is the average of the square of 
this third difference: 



N-3s 

8 S w (/)4f , . , „v 6 .... 

->U 7^7F (sm( " T/)) ' (94) 

These methods can be applied to cases in which there is dead time between 
measurements of average frequency during the sampling intervals. Suppose for 
example that the measurements consist of intervals of length r = stq during 
which an average frequency is measured, separated by dead time intervals of 
length D — t during which no measurements are available. Let the index j 
label the measurement intervals with j = 1,2, ...N. An appropriate variance 
can be defined in terms of the difference between the average frequency in the 
j th interval and that in the interval labeled by j + r: 

tfl = fyvi + T,.-Vj,.)> ( 95 ) 

where y^ s is the average frequency in the interval j of length st . Then an 
appropriate variance can be defined as 

*(^) = ((AgU 2 )- (96) 

If the measurements are sufficiently densely spaced that it is possible to pass to 
an integral, this can be shown to reduce to 

*(t,D) = ^ j\f^jl{ S m{-nfrD))\sm{-nfr)) 2 . (97) 

When D = t and r = 1 there is no real dead time and this variance reduces to 
the ordinary Allan variance. 



9 Summary and Conclusion 

In this paper a method of simulating time series for the common power-law 
noises has been developed and applied to several variances used to character- 
ize clock stability. These include overlapping and non-overlapping forms of the 
Allan variance, and the Modified Allan variance. Diagonalization of quadratic 
forms for the average variances leads to expressions for the probabilities of ob- 
serving particular values of the variance for a given sampling time r = st . 
The probabilities are expressed in terms of integrals depending on the eigen- 
values of matrices formed from squares of the second differences that are used 
to define the variance. Generally speaking, the number of eigenvalues is equal 
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to the number of terms occurring in the sum used to define averages of the 
second-difference operator, and this number gets smaller as the sampling time r 
gets larger. The probability distribution P(A) for some variance A is useful in 
estimating the ±25% confidence interval about the average variance. The eigen- 
values are usually distinct; only for white PM have eigenvalues been observed 
to occur with multiplicities other than unity. This must happen in order that 
chi-square probability distributions result. Methods for computing the prob- 
abilities have been presented in a few useful cases. Results presented in this 
paper are confined to various forms of the Allan variance. 

Other methods of simulating power law noise have been published; the 
present approach differs from that of[5] in that no causality condition is im- 
posed. The application of the methods developed here will be applied to other 
variances in future papers. 



10 Appendix 1. Evaluation of contour integrals 
for probability 

In almost all cases except for white PM, the eigenvalues are all distinct. We 
show here how then the probability function, Eq. (|43|) . can be reduced to a sum 
of real integrals. For each of the eigenvalues, we introduce the quantity 

r k = l/(2e fc ) . (98) 

The integral becomes 

r^T r a. . „iusA 

(99) 



(cj - irk) 1 / 2 



By Jordan's lemmajll). the contour of integration can be deformed into the 
upper half plane. The addition of a circular arcs at radius \uj\ = oo contributes 
nothing. Each of the square root factors in Eq. has a branch point at 

uj — ir k with a branch line extending to -Moo along the imaginary axis. We 
define the complex argument of each such factor by 

- y < arg(w - ir k ) < | . (100) 



All branch lines extend along the positive y— axis to infinity. The largest 
eigenvalues give singularities closest to the real axis. Fig. 3 illustrates the 
resulting contour. Around each branch point is a circular portion which con- 
tributes nothing because as the radius 5 of each circle approaches zero, the 
contribution to the integral approaches zero as \/S. The integral then consists 
of straight segments where uj = iy ± 5, where i5 approaches zero. Two such 
straight segments are illustrated in Fig. 3, one on each side of the branch line 
along the y— axis. Suppose the interval of interest is placed so that out of a 
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Imjoj} ReM 



Figure 3: Contour deformed to run along branch line on the y— axis. A pair of 
segments is shown with an odd number of singularities below the segments. The 
branch line required by square root factors of the form y/ui — is defined by 
the angles-either 7r/2 or — 37r/2-of the directed segments relative to the Re(w) 
axis. Contributions from pairs of segments that are above an even number of 
branch points cancel out. 

total of M eigenvalues, n of them are below the interval (three are shown in the 
figure). For the contribution to the integral on the left, y is decreasing and we 
can account for this by reversing the limits and introducing a minus sign. The 
phase factor contributed by n factors with branch points below the interval is 




(101) 



The contribution from branch points above the interval of interest is the same 
on both sides of the y— axis, and is 




(102) 



The factor in front of the integral sign in Eq. (|99| is the same on both sides of 
the y— axis and is 




(103) 



The total phase factor of this contribution, including a factor i that comes from 
setting lu = iy is thus 




(104) 
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For the integral along the segment on the positive side of the y— axis, the only 
difference is that the phase of each of the n contributions from branch points 
below the interval changes from 37r/4 to — ir/4. The phase factor for this part 
of the contour is thus 

/ \ M / \ M —n / \ n / \ n 

+ i e"T ( e M i 6 ^ 1 ) =+M e ^ Li ) ■ ( 105 ) 

If n is even, the two contributions cancel. If n = 2m + 1 is odd, then the 
contributions add up with a factor 2(— l) m . The probability is thus always real 
and consists of contributions with alternating signs, with every other interval 
left out. 

In summary, the contour integral contributions from portions of the imag- 
inary axis in the complex u> plane that have an even number of branch points 
below the interval will not contribute to the integral. For example, if there are 
four distinct eigenvalues the probability will reduce to 

P(A) - V ri r2r3r 4 f T2 e~ yA dy 



y/(y~ r i)( r 2 - y)(n - y)(ri - y) 

V(y - r i)(y - r z)(y - r 3)( r 4 - y) 



(106) 



Such results have been used to evaluate the probabilities for certain sampling 
intervals for flicker fm noise in Sect. 4. 



11 Appendix 2. Proof that Eq. (1521) generates 
positive eigenvalues 

In this Appendix we show that under all circumstances Eq. (|56[) gives rise to 
positive eigenvalues. 

Obviously if TV is a prime number Eq. (|56|) can never be satisfied (n = 0). 
Consider the case s = N — 1. Then m — MN and there are no solutions since 
m < N/2. In general both m and s may contain factors that divide M or JV. 
Suppose s = ab where at least one of the factors a, b is greater than 1, where a 
divides M and b divides N. Then let 



M = am a , N — bn b . (107) 

Then m = m a nb. m = m a rib > N/2 = biib/2 cannot happen because m < N/2. 
If m a 7ib = N/2 — brib, then the number of solutions is n = 1 and 2(s — n) — 1 = 
2(ab - 1) - 1 > 2(2 - 1) - 1 > 0. If m = m a n b < N/2 = bn b /2, then there may 
be solutions m — m a n b ,m — 2m a n b ,m — Zm a n b , ...up to bn b /2. The number 
of such solutions is 

n=L-^H = L-J (108) 

m a n b m a 
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where [^J means the largest integer less than or equal to x. The excess of 
conditions is then 

2(s-n) - 1 =2(ab- I— J) > 1, (109) 
m a 

which proves the assertion. 
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